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Abstract 

Consider a point set T> with a measure function /i : D — * M. Let A be the set of subsets of D 
induced by containment in a shape from some geometric family (e.g. axis-ahgned rectangles, 
half planes, balls, /c-oriented polygons). We say a range space {T>,A) has an e-approximation 
Pif 



max 



< £. 



We describe algorithms for deterministically constructing discrete e-app- roximations for 
continuous point sets such as distributions or terrains. Furthermore, for certain families of 
subsets A, such as those described by axis-aligned rectangles, we reduce the size of the e- 
approximations by almost a square root from 0(jt log -j) to polylogi). This is often 
the first step in transforming a continuous problem into a discrete one for which combinato- 
rial techniques can be applied. We describe applications of this result in geo-spatial analysis, 
biosurveillance, and sensor networks. 
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1 Introduction 



Representing complex objects by point sets may require less storage and may make computation on 
them faster and easier. When properties of the point set approximate those of the original object, 
then problems over continuous or piecewise-linear domains are now simple combinatorial problems 
over point sets. For instance, when studying terrains, representing the volume by the cardinality of 
a discrete point set transforms calculating the difference between two terrains in a region to just 
counting the number of points in that region. Alternatively, if the data is already a discrete point 
set, approximating it with a much smaller point set has applications in selecting sentinel nodes 
in sensor networks. This paper studies algorithms for creating small samples with guarantees in 
the form of discrepancy and e-approximations, in particular we construct e-approximations of size 
O(ipolylogi). 



e-approximations. In this paper we study point sets, which we call domains and we label as T>, 
which are either finite sets or are Lebesgue-measureable sets. For a given domain V let A be a 
set of subsets of D induced by containment in some geometric shape (such as balls or axis-aligned 
rectangles). The pair (D, A) is called a range space. We say that P is an e-approximation of (D, A) 
if 

IRnPl \Rr\V\ 



max 



\P\ |D| 



where | • | represents the cardinality of a discrete set or the Lebesgue measure for a Lebesgue- 
measurable set. A is said to shatter a discrete set X C D if each subset of X is equal to RD X 
for some R ^ A. The cardinality of the largest discrete set X that A can shatter is known as the 
VC-dimension. A classic result of Vapnik and Chervonenkis [30] states that for any range space 
{T),A) with constant VC-dimension v there exists a subset P c T> consisting ofO{-p log |) points 
that is an e-approximation for (D,7l). Furthermore, if each element of P is drawn uniformly at 
random from D such that \P\ = log then P is an e-approximation with probability at 
least 1 — 6. Thus, for a large class of range spaces random sampling produces an e-approximation 
of sizeO(4logi). 



Deterministic construction of e-approximations. There exist deterministic constructions for 

e-approximations. When D is the unit cube [0, 1]*^ there are constructions which can be interpreted 
as e-approximations of size 0( ^2d/("d+i) ) for half spaces [17] and 0{ ^2d/ld+i) log'^/*-'^^^^ ^ polylog(log 
for balls in d-dimensions [6]. Both have lower bounds of ^( ^2d/(d+i) 

) [3]. See Matousek [18] for 

more similar results or Chazelle's book [10] for applications. For a domain D, let Di^ describe the 
subsets induced by axis-parallel rectangles in d dimensions, and let describe the subsets induced 
by fc-oriented polygons (or more generally polytopes) with faces described by k predefined normal 
directions. More precisely, for /3 = . . . , l3k} C S''"^, let Qf^ describe the set of convex poly- 
topes such that each face has an outward normal ib/3j for /3j G /3. If /3 is fixed, we will use to 
denote Qp since it is the size k and not the actual set (3 that is important. When D = [0, 1]*^, then the 
range space {Uj'Jlci) has an e-approximation of size log'^^^ ^ polylog(log i)) [13]. Also, for 
all homothets (translations and uniform scalings) of any particular Q G Q^, Skriganov constructs 
an e-approximation of size 0(i log*^"^ ^ polylog(log i)). When D is a discrete point set of size 

n, e-approximations of size 0{{j log exist for bounded VC-dimension v [20], and can 
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be constructed in time 0{n ■ ^ log^ i). In this spirit, for OI2 and a discrete point set of size n, 
Suri, Toth, and Zhou [28] construct an e-approximation of size log(en) log'^(^ log(en))) in the 
context of a streaming algorithm which can be analyzed to run in time 0{n{^ log i)^). 

Our results. We answer the question, "for which ranges spaces can we construct e-approximations 
of size 0(7 polylog^)?" by describing how to deterministically construct an e-approximation of 
size polylogi) for any domain which can be decomposed into or approximated by a finite set 
of constant-size polytopes for families "Rd and Q^. In particular: 

• For a discrete point set D of cardinality n, we give an algorithm for generating an e-approximation 
for (D, Qjt) of size log^^ ^ polylog(log i)) in 0{n^ polylog^) time. This requires a 
generaUzation of the iterative point set thinning algorithm by Chazelle and Matousek [II] 
that does not rely on VC-dimension. This implies similar results for as well. 

• For any d-dimensional domain D that can be decomposed into n /c'-oriented polytopes, we 
give an algorithm for generating an e-approximation of size 0{(k+k')^\o^^ ^ polylog(log i)) 
for (©, Qfe) in time 0{{k + k')n^ polylog^). 

We are interested in terrain domains 2) defined to have a base B (which may, for instance, be a 
subset of M^) and a height function h : B ^M.. Any point (p, z) such that p e B and < z < h{p) 
(or > z > h{p) when h{p) < 0) is in the domain D of the terrain. 

• For a terrain domain D where B and h are piecewise-Iinear with n linear pieces, our result 
implies that there exists an e-approximation of size 0{k^ log^ ^ poIyIog(log ^)) for (D, Q^), 
and it can be constructed in 0(n • pr polylog^) time. 

• For a terrain domain © where S C is a rectangle with diameter d and h is smooth 

(C^ -continuous) with minimum height and largest eigenvalue of its Hessian A, we give an 
algorithm for creating an e-approximation for (D, 3^2 x I^) of size log^ ^ poIyIog(log ^)) 
in time polylogi). 

These results improve the running time for a spatial anomaly detection problem in biosurveil- 
lance [1], and can more efficiently place or choose sentinel nodes in a sensor network, addressing 
an open problem [23]. 

Roadmap. We introduce a variety of new techniques, rooted in discrepancy theory, to create e- 
approximations of size polylog^) for increasingly difficult domains. First, Section 2 discusses 
Lebesgue and combinatorial discrepancy. Section 3 generalizes and improves a classic technique 
to create an e-approximation for a discrete point set. Section 4 describes how to generate an e- 
approximation for a polygonal domain. When a domain can be decomposed into a finite, disjoint 
set of polygons, then each can be given an e-approximation and the union of all these point sets 
can be given a smaller e-approximation using the techniques in Section 3. Section 5 then handles 
domains of continuous, non-polygonal point sets by first approximating them by a disjoint set of 
polygons and then using the earlier described techniques. Section 6 shows some applications of 
these results. 
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2 Lebesgue and Combinatorial Discrepancy 



Lebesgue discrepancy. The Lebesgue discrepancy is defined for an n-point set P C [0, 1]*^ 
relative to the volume of a unit cube [0, 1]*^. ^ Given a range space ([0, 1]*^, ^1) and a point set P, the 
Lebesgue discrepancy is defined 

D{P,A) = fiViY>\D{P,R)\, where D{P,R) = n ■ \Rr\[Q,\Y\ - \Rr\ P\. 

Optimized over all n-point sets, define the Lebesgue discrepancy o/ ([0, 1]'^, A) as 

D{n,A)= inf D{P,A). 

PC[0,l]'^,\P\=n 

The study of Lebesgue discrepancy arguably began with the Van der Corput set C„ [29], which 
satisfies D{Cn,'^2) = O(logn). This was generalized to higher dimensions by Hammersley [14] 
and Halton [13] so that -D(C„, "Rd) = 0{\og'^'~^ n). However, it was shown that many lattices also 
provide 0(log n) discrepancy in the plane [18]. This is generahzed to 0{\og'^~^ n log^"''^ log n) for 
r > over "X^ [24, 25, 7]. For a more in-depth history of the progression of these results we refer 
to the notes in Matousek's book [18]. For application of these results in numerical integration see 
Niederreiter's book [21]. The results on lattices extend to homothets of any G for O(logn) 
discrepancy in the plane [24] and 0(log^^^ nlog^^^logn) discrepancy, for r > 0, in [26], 
for some constant k. A wider set of geometric families which include half planes, right triangles, 
rectangles under all rotations, circles, and predefined convex shapes produce 0(n^/^) discrepancy 
and are not as interesting from our perspective. 

Lebesgue discrepancy describes an e-approximation of ( [0, 1] , ^1) , where e = f{n) = D{n,A)/n. 
Thus we can construct an e-approximation for ([0, 1]*^, ^l) of size 5d(£, A) as defined below. (Solve 
for n in £ = D(n,A)/n).) 

^ I Oil log" i Polylog(log D) for D{n, A) = 0(log^ n) 
' \0((l/£)V(i-^)) for D{n,A) = 0{n^) 

Combinatorial discrepancy. Given a range space {X, A) where X is a finite point set and a 
coloring function x:X—*^{— 1,-1-1} we say the combinatorial discrepancy of {X,A) colored by 

disC;^(X, A) = maxdisc^(X n R) where 

disc^(X) = ^ x{x) = \{xeX : x{x) = +1}| - \{x e X : x(x) = -1}| . 
xex 

Taking this over all colorings and all point sets of size n we say 

disc(n,yi) = max min disCy(X,yi). 

{X:\X\=n}x:X^{-l,+l} 

Results about combinatorial discrepancy are usually proved using the partial coloring method [5] 
or the Beck-Fiala theorem [9]. The partial coloring method usually yields lower discrepancy by 



'Although not coimnon in the literature, this definition can replace [0, 1]'' with an hyper-rectangle [0, wi] x [0, W2] x 
X [0,Wd]. 
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some logarithmic factors, but is nonconstructive. Alternatively, the Beck-Fiala theorem actually 
constructs a low discrepancy coloring, but with a sUghtly weaker bound. The Beck-Fiala theorem 
states that for a family of ranges A and a point set X such that max^jgx |{74Gyi : x ^ A}\ < t, 
disc(X, A) <2t—l.So the discrepancy is only a constant factor larger than the largest number of 
sets any point is in. 

Srinivasan [27] shows that disc(n, 3^2) = O(log^'^n), using the partial coloring method. An 
earUer result of Beck [4] showed disc(n, = 0(log'^ n) using the Beck-Fiala theorem [9]. The 
construction in this approach reduces to 0(n) Gaussian eUminations on a matrix of constraints that 
is 0(n) X 0{n). Each Gaussian elimination step requires 0{n?) time. Thus the coloring x i^i 
the construction for disc(n, 3^2) = 0{\og* n) can be found in 0(n*^) time.We now generalize this 
result. 

Lemma 2.1. disc(n, Q^) = 0(log^'^ n) for points in and the coloring that generates this dis- 
crepancy can be constructed in 0{n^) time, for k constant. 

The proof combines techniques from Beck [4] and Matousek [19]. 

Proof. Given a class Qfc, each potential face is defined by a normal vector from {f3\, . . . , f5k}- For 
G [1, k] project all points along (5j. Let a canonical interval be of the form for integers 

q G [l,logn] and t G [0,2'^). For each direction j3j choose a value q G [l,logn] creating 2^ 
canonical intervals induced by the ordering along f3j . Let the intersection of any k of these canonical 
intervals along a fixed f3j be a canonical subset. Since there are log n choices for the values of q for 
each of the k directions, it follows that each point is in at most (logn)'^ canonical subsets. Using 
the Beck-Fiala theorem, we can create a coloring for X so that no canonical subset has discrepancy 
more than 0{\.og^ n). 

Each range G Qfc is formed by at most ©(log*^ n) canonical subsets. For each ordering by 
f3i, the interval in this ordering induced by R can be described by O(logra) canonical intervals. 
Thus the entire range R can be decomposed into 0(log'^ n) canonical subsets, each with at most 
0(log^ n) discrepancy. 

Applying the Beck-Fiala construction of size n, this coloring requires 0{n*) time to construct. 

□ 

Corollary 2.1. disc{n, = 0(log^'^ n) and the coloring that generates this discrepancy can be 
constructed in O(n^) time, for d constant. 

A better nonconstructive bound exists due to Matousek [19], using the partial coloring method. 
For polygons in disc(n, Q^) = 0{k log^'^ nY^log(A; + log n)), and for polytopes in M.'^ disc(n, Qk) = 
Q^j^i.5ld/2\ log'^+^Z^ n-yiog(fe~+Togn)). For more results on discrepancy see Beck and Chen's 



Similar to Lebesgue discrepancy, the set P = {p G X \ x{p) = +1} generated from the 
coloring x for combinatorial discrepancy disc(n. A) describes an e-approximation of {X, A) where 
e = f{n) = disc(n,yi)/n. Thus, given this value of s, we can say that P is an e-approximation for 
{X, A) of size 



The next section will describe how to iteratively apply this process efficiently to achieve these 
bounds for any value of e. 



book [8]. 




(2.2) 
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3 Deterministic Construction of ^-approximations for Discrete 
Point Sets 



We generalize the framework of Chazelle and Matousek [11] describing an algorithm for creating 

an e-approximation of a range space {X,A). Consider any range space {X,A), with \X\ = n, 
for which there is an algorithm to generate a coloring x that yields the combinatorial discrepancy 
disc^(X,A) and can be constructed in time 0{n^ ■ l{n)) where l{n) = o(n). For simplicity, we 
refer to the combinatorial discrepancy we can construct disc^(X, A) as disc(n, A) to emphasize the 
size of the domain, and we use equation (2.2) to describe g{e, A), the size of the e-approximation it 
corresponds to. The values disc(n, A), w, and l{n) are dependent on the family A (e.g. see Lemma 
2.1), but not necessarily its VC-dimension as in [11]. As used above, let /(n) = disc(n,>l)/n be 
the value of e in the e-approximation generated by a single coloring of a set of size n — the relative 
error. We require that, /(2n) < (1 — (5)/(n), for constant < S < 1; thus it is a geometrically 
decreasing function. 

The algorithm will compress a set X of size n to a set P of size 0{g{e,A)) such that P is an 
e-approximation of {X, A) by recursively creating a low discrepancy coloring. We note that an 
£-approximation of an e'-approximation is an (e -|- £')-approximation of the original set. 

We start by dividing X into sets of size 0{g{e,A)),^ here e is a parameter. The algorithm 
proceeds in two stages. The first stage alternates between merging pairs of sets and halving sets by 
discarding points colored xip) = —1 by the combinatorial discrepancy method described above. 
The exception is after every w + 2 halving steps, we then skip one halving step. The second stage 
takes the one remaining set and repeatedly halves it until the error f{\P\) incurred in the remaining 
set P exceeds 2+26- results in a set of size 0{g{e, A)). 



Algorithm 3.1 Creates an e-approximation for {X,A) of size 0{g{e,A)). 

1: Divide X into sets {Xq, Xi , X2, . . .} each of size 4(w -|- 2)g{e, A). ^ 

2: repeat {Stage 1} 

3: for w + 2 steps do {or stop if only one set is left} 

4: Merge: Pair sets arbitrarily (i.e. Xi and Xj) and merge them into a single set (i.e. 
Xi := Xi U Xj). 

5: Halve: Halve each set Xi using the coloring x from disc{Xi,A) (i.e. Xi = {x e Xi \ 

X{x) = +1}). 

6: end for 

7: Merge: Pair sets arbitrarily and merge each pair into a single set. 

8: until only one set, P, is left 

9: repeat {Stage 2} 

10: Halve: Halve P using the coloring x from disc(P, A). 

11: until /(|P|) > e/(2 + 26) 



Theorem 3.1. For a finite range space {X,A) with \X\ = n and an algorithm to construct a 
coloring x '■ X — > {—1,-1-1} such that 

^If the sets do not divide equally, artificially increase the size of the sets when necessary. These points can be removed 
later. 
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• the set {x £ X : xi^) = +1} '■^ a-approximation of{X,A) of size g{a,A) with 
a = disc-^{X,A)/n (see equation (2.2)). 

• X c«n be constructed in 0{n^ ■ l{n)) time where l{n) = o{n). 

then Algorithm 3.1 constructs an e-approximation for (X, ^l) of size 0{g{e, A)) in time 0{w'^~^n ■ 
g{e,Ar-^-l{g{e,A))+g{e,A)). 

Proof. Let 2^ = 4:{w + '2)g{£, A), for an integer j, be the size of each set in the initial dividing stage 
(adjusting by a constant if 5 < \). Each round of Stage 1 performs w; + 3 MERGE steps and w + 2 
Halve steps on sets of the same size and each subsequent round deals with sets twice as large. 
The union of all the sets is an o-approximation of (X, A) (to start a = 0) and a only increases in 
the Halve steps. The i\h round increases a by /(2-'""^+*) per Halve step. Since /(n) decrease 
geometrically as n increases, the size of a at the end of the first stage is asymptotically bounded 
by the increase in the first round. Hence, after Stage 1 a < 2{w + 2)/(4(it; + 2)g{e,A)) < |. 
Stage 2 culminates the step before fi\P\) > Thus the final Halve step creates an 

approximation and the entire second stage creates an | -approximation, hence overall Algorithm 3.1 
creates an e-approximation. The relative error caused by each Halve step in stage 2 is equivalent 
to a Halve step in a single round of stage 1. 

The running time is also dominated by Stage 1. Each Halve step of a set of size 2^ takes 
0{{2^)'^l{2^)) time and runs on n/2^ sets. In between each HALVE step within a round, the number 
of sets is divided by two, so the running time is asymptotically dominated by the first Halve step 
of each round. The next round has sets of size 2-'+^, but only n/2^^'^^'^ of them, so the runtime 
is at most ^ that of the first Halve step. Thus the running time of a round is less than half of that 
of the previous one. Since 2^ = 0{wg{£,A)) the running time of the HALVE step, and hence the 
first stage is bounded by 0{n ■ (w ■ g{s,A))'^~^ ■ l{g{£,A)) + g{£,A)). Each Halve step in the 
second stage corresponds to a single Halve step per round in the first stage, and does not affect the 
asymptotics. □ 

We can invoke Theorem 3.1 along with Lemma 2.1 and Corollary 2.1 to compute x in 0{n'^) 
time (notice that w = i and /(•) is constant), so g{£,Qk) = 0(i log^'^ ^ poly log (log ^)) and 
g{e,'Jld) =0(7 log^^ \ poly log (log i)). We obtain the following important corollaries. 

Corollary 3.1. For a set of size n and over the ranges Qk an e-approximation of size log^^ ^ polylog(log ^)) 
can be constructed in time 0{n^ polylog^). 

Corollary 3.2. For a set of size n and over the ranges "R^ e-approximation of size O ( ^ log^'' ^ polylog (log ^ ) ) 
can be constructed in time 0{n^ polylog i). 

Weighted case. These results can be extended to the case where each point x e X is given a 
weight n{x). Now an e-approximation is a set P C X and a weighting /j, : X ^M. such that 



max 



n{P n R) ii{x n R) 



m(p) KX) 



where /i(P) = J2peP l^iP)- ^he weights on P may differ from those on X. A result from Matousek 
[16], invoking the unweighted algorithm several times at a geometrically decreasing cost, creates 
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a weighted e-approximation of the same asymptotic size and with the same asymptotic runtime 
as for an unweighted algorithm. This extension is important when we combine e-approximations 
representing regions of different total measure. For this case we weight each point relative to the 
measure it represents. 

4 Sampling from Polygonal Domains 

We will prove a general theorem for deterministically constructing small £-approximations for 
polygonal domains which will have direct consequences on polygonal terrains. A key observation of 
Matousek [16] is that the union of e-approximations of disjoint domains forms an e-approximation 
of the union of the domains. Thus for any geometric domain D we first divide it into pieces for which 
we can create e-approximations. Then we merge all of these point sets into an e-approximation for 
the entire domain. Finally, we use Theorem 3.1 to reduce the sample size. 

Instead of restricting ourselves to domains which we can divide into cubes of the form [0, 1]'', 
thus allowing the use of Lebesgue discrepancy results, we first expand on a result about lattices and 
polygons. 

Lattices and polygons. For a; G M, let \x[ represent the fractional part of x, and for a G M'^^^ 
let a = (ofi , . . . , ccrf- 1 ) . Now given a and m let Pa,m = {po > ■ ■ ■ > Pm- 1 } be a set of m lattice points 
in [0, l]'^ defined pi = {^^,\aii[, . . . , Ja^^iil.). Pa,m is irrational with respect to any polytope in 
Qg if for all /3j € (3, for all j < d, and for sM h < d — \, the fraction [^ij/ah is irrational. (Note 
that /3ij represents the jth element of the vector f3i.) Lattices with a irrational (relative to the face 
normals) generate low discrepancy sets. 

Theorem 4.1. Let Q G Q^/ be a fixed convex polytope. Let /3, j3' C S*^"^ be sets of k and k' direc- 
tions, respectively. There is an e-approximation of{Q, Qjj) of size 0{{k-\-k')^ log'^"^ ^ polylog(log ^)). 

This e-approximation is realized by a set of lattice points Pa,m H Q such that Pa,m is irrational 
with respect to any polytope in Q/3u/3'- 

Proof. Consider polytope tQh and lattice Pa,m, where the uniform scaling factor t is treated as an 
asymptotic quantity. Skriganov's Theorem 6.1 in [26] claims 



for r > 0, as long as Pa,m is irrational with respect to the normal of the face / of and infinite 
otherwise, where 6* G (0, 1) and p can be arbitrarily large. Note that this is a simphfied form yielded 
by invoking Theorem 3.2 and Theorem 4.5 from [26]. By setting p^ = t'^~^. 



max 




where 



Sf{Pa,m,p) = 0(l0g'^-VW+"l0gp) 



maxD{Pa,m, tQh + v)= 0{h\og'^-Hlog^+^ logt). 



(4.1) 
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Now by noting that as t grows, the number of lattice points in tQh grows by a factor of t , and we can 
sett = n^/'^ so(4.1)impUes thatD(Pa,TO,i(5/i) = 0(/ilog°'"^ nlog^+^logn) for \Pa,m\ = m = n 
and tQh C [0, 1]'^. 

The discrepancy is a sum over the set of h terms, one for each face /, each of which is small as 
long as Pa,m is irrational with respect to /'s normal Hence this lattice gives low discrepancy 
for any polytope in the analogous family such that Pa,m is irrational with respect to Q^. Finally 
we realize that any subset QCiQ^for Q e Qp' and Qk G Q/3 is a polytope defined by normals from 

U /? and we then refer to 5£)(£, Q/3u/3') i" (^-l) to bound the size of the e-approximation from the 
given Lebesgue discrepancy. □ 

Remark. Skriganov's result [26] is proved under the whole space model where the lattice is 
infinite (tQh is not confined to [0, l]'^), and the relevant error is the difference between the measure 
of tQh versus the cardinality \tQh n Pa,m\y where each p e Pa,Tn represents 1 unit of measure. 
Skriganov's main results in this model is summarized in equation (4.1) and only pertains to a fixed 
polytope Qh instead of, more generally, a family of polytopes Qp, as shown in Theorem 4.1. 

Samples for polygonal terrains. Combining the above results and weighted extension of The- 
orem 3.1 impUes the following results. 

Theorem 4.2. We can create a weighted e-approximation of size 0{{k+k')^-\og^'' ^ polylog(log i)) 
of{D, Qfe) in time 0{{k + k')n-^ polylogi)/or any d-dimensional domain D which can be decom- 
posed into n d-dimensional convex k' -oriented polytopes. 

Proof. We divide the domain into n /c'-oriented polytopes and then approximate each polytope Qk' 
with a point set Pa.m H Qk' using Theorem 4.1. We observe that the union of these point sets is a 
weighted e-approximation of (D, Q^), but is quite large. Using the weighted extension of Theorem 
3. 1 we can reduce the point sets to the size and in the time stated. □ 

This has applications to terrain domains V defined with a piecewise-linear base B and height 
function h : B ^ M.. We decompose the terrain so that each linear piece of h describes one 
3-dimensional polytope, then apply Theorem 4.2 to get the following result. 

Corollary 4.1. For terrain domain D with piecewise-linear base B and height function h : B ^M. 
with n linear pieces, we construct a weighted e-approximation of{T>, Q^) of size 0(k^ log^ ^ polylog(lo 
in time 0{kn-^ polylog^). 

5 Sampling from Smooth Terrains 

We can create an e-approximation for a smooth domain (one which cannot be decomposed into 
polytopes) in a three stage process. The first stage approximates any domain with a set of polytopes. 
The second approximates each polytope with a point set. The third merges all point sets and uses 
Theorem 3.1 to reduce their size. 

This section mainly focuses on the first stage, however, we also offer an improvement for the 
second stage in a relevant special case. More formally, we can approximate a non-polygonal domain 
D with a set of disjoint polygons P such that P has properties of an e-approximation. 
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Lemma 5.1. If \T> \ P\ < ^\D\ and P Q D then max 



\Rr\P\ \Rr\D\ 



\P\ 



< e. 



Proof. No range R ^ A can have 

I'i^li 



l-KnP| 



©I 



|i?nD| 

m 



> £ because if |D| > |Pj (w.l.o.g.), then 

|D| 



i?nXi| -|^|i?nP| < e\T)\ and |i2nP|^ - |i?nD| < e|T)|. The first part follows from jf| > 1 



and is loose by a factor of 2. For the second part we can argue 



\RnPY--p--\Rr\T)\ < 



\Rnp\ 



1 - - 

^ 2 



1 - - 

^ 2 



iRnDl < slRnvi < em. 



□ 



For terrain domains D defined with a base B and a height function /t : S — >^ R, if S is polygonal 
we can decompose it into polygonal pieces, otherwise we can approximate it with constant-size 
polygonal pieces according to Lemma 5.1. Then, similarly, if h is polygonal we can approximate 
the components invoking Corollary 4. 1 ; however, if it is smooth, then we can approximate each 
piece according to Lemma 5.1. 

Section 5.1 improves on Theorem 4.1 for the second stage and gives a more efficient way to create 
an ^-approximation for (D, 3?^ x ]R) of a terrain when i? is a rectangle and h is linear. Ranges from 
the family 31^ x M are generalized hyper-cylinders in d -|~ 1 dimensions where the first d dimensions 
are described by an axis-parallel rectangle and the {d + l)st dimension is unbounded. Section 5.2 
focuses on the first stage and uses this improvement as a base case in a recursive algorithm (akin to 
a fair spUt tree) for creating an e-approximation for (D, 31^ x R) when B is rectangular and h is 
smooth. 



5.1 Stretching the Van der Corput Set 

The Van der Corput set [29] is a point set P„ = {po, . . . ,Pn-i} in the unit square defined for 
Pi = {^,b{i)) where b{i) is the bit reversal sequence. For simplicity we assume n is a power of 2. 
The function b{i) writes i in binary, then reverses the ordering of the bits, and places them after the 
decimal point to create a number in [0, 1). For instance forn = 16, i = 13 = 1101 in binary and 
6(13) = 0.1011 = ii. Formally, if i = ai2' then 5(i) = E'T 

Halton [13] showed that the Van der Corput set Pn satisfies D(P„,3?2) = O(logn). We can 
extend this to approximate any rectangular domain. For a rectangle [0, w] x [0, 1] (w.l.o.g.) we can 
use the set Pn,w,i where pi = {w ■ ^,1 ■ b{i)) and a version of the Lebesgue discrepancy over a 
stretched domain is still O(logra). 

We can stretch the Van der Corput set to approximate a rectangle r = [0, w] x [0, 1] with a 
weighting by an always positive linear height function h{x, y) = ax + (5y + 7. Let A('u;, a, 7, i) 
be defined such that the following condition is satisfied 

A(?i),a,7,i) ^ 

(ax + j)dx = — (ax + ^)dx. 
n Jo 

Note that we can solve for A explicitly and because h is linear it can simultaneously be defined 
for the X and y direction. Now define the stretched Van der Corput set Sn,w,l,h = {so, ■ ■ ■ , ■Sn} for 
Si = {A{w, a, 7, i),A{l, p, 7, b(i) ■ n)). 
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Theorem 5.1. For the stretched Van der Corput set Sn,w,i,h, D{Sn.w.i.h, '^2) = O(logn) over the 
domain [0, w] x [0, 1] with h : [0, w\ x [0, 1] M"*" a linear weighting function. 

The proof follows the proof inMatousek [18] for proving logarithmic discrepancy for the standard 
Van der Corput set in the unit square. 



A(f,/3,7,fc) A(;,/3,7,A:+ 
29 ' 21 



j for integers q G [1, n] and 



Proof. Let a canonical interval be of the form 

k G [0,2'?). Let any rectangle r = [0, a) x / where / is canonical and a G (0, 1] be called a 
canonical rectangle. 

Claim 5.1. For any canonical rectangle r, D{Sn,w,l,hi f) < 1- 

Proof. Like in the Van der Corput set, every subinterval of r such that h{r) = ^ has exactly 1 point. 

Let I = [^Ml^, ^M^y Thus each rectangle r,- = [A(Z, ^, 7, A(/, 7, ^^)) x 

/ contains a single point from Sn,w,i,h and = ^, where h{r) = h{p)dp. 

So the only part which generates any discrepancy is the canonical rectangle rj which contains the 
segment a x I. But since \Sn^w,l,h n r^- fl r| < 1 and h{rj n r) < i, the claim is proved. □ 

Let Crf be the family of ranges consisting of d-dimensional rectangles with the lower left comer 
at the origin. Let C(x,j/) € C2 be the comer rectangle with upper right comer at {x, y). 

Claim 5.2. Any comer rectangle Cf^^^y-^ can be expressed as the disjoint union of at most 0(log n) 
canonical rectangles plus a rectangle M with \D{Sn,w,l,h-i ^)\ < 1- 

Proof. Starting with the largest canonical rectangle rg = [0, a) x I within C^^^y) such that / = 
A(^,A7.o) ^ A(;,^^,7,i) ^ smallest value possible of q, keep including the next largest disjoint 

canonical rectangle within C(^x,y)- Each consecutive one must increase q by at least 1. Thus there 
can be at most 0(log n) of these. 

The left over rectangle M = [nix , x] x [my , y] , must be small enough such that h{p, q)dqdp < 

^, thus it can contain at most 1 point and D{Sn,w,l,hj < 1. □ 

It follows from Claim 5.1 and Claim 5.2 that disc(S', 62) = 0(log n). We conclude by using the 
classic result [18] that D{S, 62) < D{S, ^2) < 41^(5, 62) for any point set S. □ 

This improves on the discrepancy for this problem attained by using Theorem 4. 1 by a factor of 
logi. 

Corollary 5.1. A stretched Van der Corput set Sn^w.lh ff^^^ <^fi e-approximation of (D,3?2) of 
size n = 0(i log ^ polylog(log \))for T) defined by a rectangle [0, w\ x [0, l\ with a linear height 
function h. 

Remark. This extends to higher dimensions. A stretched b-ary Van der Corput set [18] forms 
an e-approximation of (D, of size 0(i log°'^^ ^ polylog(log ^)) for D defined by x^^^[0, Wi] 
with a linear height function. Details are omitted. 
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5.2 Approximating Smooth Terrains 

Given a terrain domain D where B C is rectangular and h : B ^ is a C^-continuous height 
function we can construct an e-approximation based on a parameter e and properties dxi, and 
Ad. Let = min^gs h{p). Let d-j) = maxp^^g^ ||p — be the diameter of D. Let Ad be the 



largest eigenvalue of where = 



dx'-' dxdy 
d'^h d^h 



is the Hessian of h. 



dydx dy^ 

We first create a set of Unear functions to approximate h with a recursive algorithm. If the entire 

domain cannot be approximated with a single Unear function, then we split the domain by its longest 
direction (either x or y direction) evenly. This decreases d-ri by a factor of 1/ \/2 each time. We 
recur on each subset domain. 

Lemma 5.2. For a domain T> with rectangular base B cM? with a -continuous height function 
h : B ^ M. we can approximate h with 0{ ) linear pieces so that for all p ^ B h^{p) < 
h{p) < heip) + e. 

Proof. First we appeal to Lemma 4.2 from Agarwal et. al [2] which says that the error of a first 
order linear approximation at a distance d is bounded by A^c^^- Thus we take the tangent at the 
point in the middle of the range and this Unear (first order) approximation has error bounded by 
AD(dD/2)^ = ADdD/4. The height of the Unear approximation is lowered by X-rid'^/A from the 
tangent point to ensure it describes a subset of D. Thus, as long as the upper bound on the error 

Adc?|)/2 is less than z^e then the lemma holds. The ratio is halved every time the domain is 
split until it is less than 1. Thus it has 0(^^^) base cases. □ 

After running this decomposition scheme so that each linear piece L has error e/2, we invoke 
Corollary 5.1 to create an (e/2)-approximation point set of size log ^ polylog(log i)) for each 
{L,^2 X IR)- The union creates a weighted e-approximation of {T>,Jl2 x M), but it is quite large. 
We can then reduce the size according to Corollary 3.2 to achieve the following result. 

We can improve further upon this approach using a stretched version of the Van der Corput Set 
and dependent on specific properties of the terrain. Consider the case where 5 is a rectangle with 
diameter d^ and h is continuous with minimum value z^ and where the largest eigenvalue of 
its Hessian is Ad- For such a terrain D, interesting ranges x M are generalized cylinders where 
the first 2 dimensions are an axis-parallel rectangle and the third dimension is unbounded. We can 
state the following result (proved in the full version). 

Theorem 5.2. For a domain T> with rectangular base i? C and with a -continuous height 
function h : B ^Wwe can deterministically create a weighted e-approximation of {D,^2 x M) of 

size O ^ (7 log"^ I poly log (log e))^ • reduce the size to log*^ ^ polylog(log j)) in 

rimeo((^)^polylogi). 

This generalizes in a straightforward way for B G M''. Similar results are possible when B is not 
rectangular or when B is not even piecewise-linear. The techniques of Section 4 are necessary if 
Qk is used instead of and are slower by a factor O(^). 
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6 Applications 



Creating smaller e-approximations improves several existing algorithms. 

6.1 Terrain Analysis 

After creating an e-approximation of a terrain we are able to approximately answer questions about 
the integral over certain ranges. For instance, a terrain can model the height of a forest. A foresting 
company may deem a region ready to harvest if the average tree height is above some threshold. 
Computing the integral on the e-approximation will be much faster than integrating over the full 
terrain model. 

More interesting analysis can be done by comparing two terrains. These can represent the forest 
height and the ground height or the elevation of sand dunes at two snapshots or the distribution 
of a population and a distribution of a trait of that population. Let Ti and T2 be two terrains 
defined by piecewise-linear height functions hi and h2, respectively, over a subset of M^. The 
height h = hi — h2 may be negative in some situations. This can be handled by dividing it into 
two disjoint terrains, where one is the positive parts of h and the other is the negative parts. Each 
triangle can be spUt by the ^ = plane at most once, so this does not asymptotically change the 
number of piecewise-linear sections. 

Once an e-approximation has been created for the positive and negative terrain, the algorithms 
of Agarwal et. al. [2] can be used to find the rectangle with the largest positive integral. For n 
points this takes 0(n^ logn) time. The same can be done for finding the rectangular range with 
the most negative integral. The range returned indicates the region of largest difference between the 
two terrains. The runtime is dominated by the time to create the ^-approximation in Corollary 4.1. 

6.2 Biosurveillance 

Given two points set representing measured data M and representing baseUne data B, anomaly 
detection algorithms find the region where M is most different from B. The measure of difference 
and limits on which regions to search can vary significantly [15, 1, 22]. One well-formed and 
statistically justified definition of the problem defines the region R from a class of regions A that 
maximizes a discrepancy function based on the notion of spatial scan statistics [15, 2]. Where 
mji = |i? n Af|/|M| and bn = \Rn B\/\B\ represent the percentage of the baseline and measured 
distributions in a range R, respectively, then the Poisson scan statistic can be interpreted as the 
Poisson discrepancy function dp{mR, bn) = uirIii ^ + (1 — mR) In . This has important 
appUcations in biosurveillance [15] where B is treated as a population and M is a subset which 
has a disease (or other condition) and the goal is to detect possible regions of outbreaks of the 
disease as opposed to random occurrences. We say a linear discrepancy function is of the form 
di(mR, hu) = aniR + f3bR + 7 for constants a, (5, and 7. The Poisson discrepancy function can 
be approximated up to an additive e factor with log^ n) linear discrepancy functions [2]. The 
range i? G 3?2 which maximizes a linear discrepancy function can be found in 0{v? log n) time and 
the G 3?2 which maximizes any discrepancy can be found in O(n^) time where + |M| = n. 

Agarwal et. al. [1] note that a random sample of size log ^) will create an e-approximation 
with probabiUty \ — 5. This can be improved using Corollary 3.2. We then conclude: 
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Theorem 6.1. Let \M[JB\ = n. A range i? G 3^2 such that \dp[mii, — maxrg3j2 dp(mr, br)\ < 
e can be deterministically found in 0{nj^ poly log (log ^) + ^ poly log (log i)) time. 

A range R E:^2 such that \dp{mji, bp) — maxj-gj^j dp{mr, br)\ < s + S can be deterministically 
found in 0{n^ polylog(log + polylog(log \)) time. 

This can be generalized to when M and B are terrain domains. This case arises, for example, 
when each point is replaced with a probability distribution. 



Generating Terrains with Kernel Functions. A drawback of the above approach to finding 
maximum discrepancy rectangles is that it places the boundaries of rectangles in some arbitrary 
place between samples. This stems from the representation of each sample as a fixed point. In 
reality its location is probably given with some error, so a more appropriate model would be to 
replace each point with a kernel density function. Probably, the most logical kernel function would 
be a Gaussian kernel, however, this is continuous and its tails extend to infinity. The base domain 
can be bounded to some polygon B so that the integral under the kernel function outside of B is 
less than e times the entire integral and then Theorem 5.2 can be applied. (See Appendix B for 
details.) Alternatively, we can replace each point with a constant complexity polygonal kernel, like 
a pyramid. Now we can ask questions about spatial scan statistics for a measured Tm and a baseUne 
Tb terrain. 

For simple ranges such as (axis-parallel half spaces) and S|| (axis-parallel slabs) finding the 
maximal discrepancy ranges on terrains reduces to finding maximal discrepancy intervals on points 
sets in 

Theorem 6.2. For a terrain T defined by a piecewise-linear height function h with n vertices 

arg max / h(p) dp and arg max / h(p) dp 

can be found in 0{n) time. 

Proof, (sketch) Project all terrains onto the axis perpendicular to halfplanes (or slabs). Integrate 
between points where the projected terrain crosses 0. Treat these intervals as weighted points and 
use techniques from Agarwal et. al. [2]. The full proof is given in Appendix A. □ 

However for 3^2, this becomes considerably more difficult. Under a certain model of computation 
where a set of 4 quadratic equations of 4 variables can be solved in constant time, the maximal 
discrepancy rectangle can be found in O(n^) time. However, such a set of equations would require 
a numerical solver, and would thus be solved approximately. But using Theorem 6.1 we can answer 
the question within en 'mO(n^ poly log ^ -|- ^ polylogi) time for a terrain with 0{n) vertices. 

Alternatively, we can create an e-approximation for a single kernel, and then replace each point in 
M and B with that e-approximation. Appendix B describes, for a Gaussian function if, how to cre- 
ate an e-approximation for 3I2) of size polylog(log i)) in time polylog(log i)). For 
standard kernels, such as Guassians, we can assume such e-approximations may be precomputed. 
We can then apply Corollary 3.1, before resorting to techniques from Agarwal et al. [2]. 

Theorem 6.3. Let \M\JB\ = n, where M and B are point sets describing the centers of Gaussian 

. r, m 1 fcp Mix) . r c R Mix) 

kernels with nxed variance, tor a range K G Jio, let uir = r"" \ and let or = ^ \ - 

A range i? G OI2 such that \dp{mR, bn) — max^gj^g dp(mr, br)\ < s can be deterministically 
found in 0{n^ polylog(log ^)) time. 
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6.3 Cuts in Sensor Networks 

Sensor networks geometrically can be thought of as a set of n points in a domain D. These points 
(or nodes) need to communicate information about their environment, but have limited power. Shri- 
vastava et. al. [23] investigates the detection of large disruptions to the domain that affect at least 
en nodes. They want to detect these significant events but with few false positives. In particular, 
they do not want to report an event unless it affects at least |n nodes. 
We say P C D is an e-sentinel of (D, A) if for all G A 

• if |i?nD| > £|D| then |i?nP| >£||P|, and 

• if |PnP| > e||P| then |i2n!D| > 

Shrivastava et. al. [23] construct e-sentinels for half spaces of size 0{j) and in expected time 
0(| log n). They note that an e/4-approximation can be used as an e-sentinel, but that the standard 
upper bound for e-approximations [30] requires roughly log ^) points which is often imprac- 
tical. They pose the question: For what other classes of ranges can an e-sentinel be found of size 
O(ipolylogi).^ 

FoUowup work by Gandhi et. al. [12] construct e-sentinels for any A with bounded VC- 

dimension v (such as disks or ellipses) of size log ^) and in time 0{n^ log^ ^). 

As an alternative to this approach, by invoking Corollary 3.1 we show that we can construct a 
small e-sentinel for Q^. 

Theorem 6.4. For a discrete point set D of size n, we can compute e-sentinels for (D, Q^) of size 
0(1 log2fc i polylog(log i)) in time 0{njj polylog{\og \)). 

In fact, if we can choose where we place our nodes we can create an e-sentinel of size 0(i polylog^) 
to monitor some domain D. We can invoke Theorem 4.1 or Theorem 5.2, depending on the nature 
of D. 

Additionally, by modifying the techniques of this paper, we can create 0(ne/ log^'^ ^) dis- 
joint sets of e-sentinels. At every Halve step of Algorithm 3.1 we make a choice of which 
points to discard. By branching off with the other set into a disjoint e-approximation, we can 
place each point into a disjoint e-sentinel of size 0(Mog^'^ polylog(log ^)). Since the Halve 
step now needs to be called 0(rie/ log^*^ times on each of the 0(log(ne)) levels, this takes 
0{n^ log(ne) polylog(log i)) time. 

Theorem 6.5. For a discrete point set D of size n, we can create 0{ne/ \o^^ ^) disjoint sets of 
e-sentinels in 0{n-^ log(rae) polylog(log ^)) total time. 

The advantage of this approach is that the nodes can alternate which sensors are activated, thus 
conserving power. If instead a single node is used in multiple e-sentinels it will more quickly use up 
its battery supply, and when its batter runs out, the e-sentinels using that node can no longer make 
the appropriate guarantees. 
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A Combinatorial Algorithms on Terrains 



A.l Half spaces, intervals, and slabs 

Let : M — >^ M be a piecewise-linear height function over a one-dimensional domain with a possibly 
negative range. Each range in 3-C| | is defined by a single point in the domain. 

Lemma A.l. For continuous /i : M — R the 



arg max / h{p) dp 



is defined by an endpoint r such that h{r) = 0. 

Proof. If the end point r moved so the size of R is increased and h{r) > then the integral would 
increase, so h{r) must be non positive. If the end point r is moved so the size of R is decreased and 
h{r) < then integral would also increase, so h{r) must be non negative. □ 

This proof extends trivially to axis-parallel slabs §|| (which can be thought of as intervals) as well. 
Lemma A.2. For continuous /i : M — ^ R, the 



arg max / h(p) dp 



is defined by two endpoints ri and such that h{ri) = and h{rr) = 0. 

Let h have n vertices. For both | and S 1 1 , the optimal range can be found in O (n) time. For JC| | , 
simply sweep the space from left to right keeping track of the integral of the height function. When 
the height function has a point r such that h{r) = 0, compare the integral versus the maximum so 
far. 

For §1 1 , we reduce this to a one-dimensional point set problem. First sweep the space and calculate 
the integral in between every consecutive pair of points ri and r2 such that h{ri) = = h{r2) and 
there is no point r-^ such that h{r^) = and ri < < r2. Treat each of these intervals as a 
point with weight set according to its value. Now run the algorithm from Agarwal et al. [2] for 
Unear discrepancy of red and blue points where the positive intervals have a red weight equal to the 
integral and the negative intervals have a blue weight equal to the negative of the integral. 

Theorem A.l. For continuous, piecewise-linear ^ : R — ^ M with n vertices 

arg max / h(p) dp and arg max / h(p) dp 
^Re^Jr ReSuJR 

can be calculated in 0{n) time. 

This result extends trivially to a higher dimensional domains as long as the famihes of ranges are 
no more compUcated. 

Theorem A.2. [Theorem 6.2] For continuous, peicewise-linear h with n vertices, 

rg max / h(p) dp and arg max / h(p) dp 
^Re:KjR ReS^JR 



are 



can be calculated in 0{n) time. 
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Proof. The sweep step of the algorithms described above are performed in the same way, only now 
the integral of up to 0(n) cubic functions must be calculated. However, this integral can be stored 
implicitly as a single linear function and can be updated in constant time every time a new vertex is 



Finally, we not a slightly surprising theorem about the difference between two terrains. 

Theorem A.3. Let M : M*^ — > M and B : ^ 'M.be piecewise-linear functions with n and m 
vertices respectively. 



can be calculated in 0{n + m) time. 

Naively, this could be calculated in 0{nm) time by counting the vertices on the terrain h{p) = 
M{p) — B{p). But we can do better. 

Proof. Although there are more than n + m vertices in h{p) = M{p) — B (p), the equations describ- 
ing the height functions only change when a vertex of one of the original functions is encountered. 
Thus there are only 0{n + m) linear functions which might cross 0. These can be calculated by 
projecting M and B independently to the axis of or S|| and then taking their sum between each 
consecutive vertex of either function. □ 

A. 2 Rectangles 

Although the work by Agarwal et al. [2] extends the one-dimensional case for point sets to a 
0(n^ log n) algorithm for rectangles, when the data is given as picewise-linear terrains the direct 
extension does not go through. However, a simple 0{n'^) time algorithm, under a certain model, 
does work. Following algorithm Exact from Agarwal et al. [1], we make four nested sweeps over 
the data. The first two bound the x coordinates and the second two bound the y coordinates. The 
inner most sweep keeps a running total of the integral in the range. However, unlike Exact each 
sweep does not give an exact bound for each coordinate, rather it just restricts its position between 
two vertices. The optimal position is dependent on all four positions, and needs to be determined 
by solving a system of four quadratic equations. This system seems to in general have no closed 
form solution (see the next subsection) and needs to be done via a numerical solver. However, these 
equations can be updated in constant time in between states of each sweep, so under the model that 
the numerical solver takes 0(1) time, this algorithm runs in O(n^) time. 

For the full correctness of the algorithm, there is actually one more step required. Given that each 
side of the rectangle is bounded between two vertices, the set of four equations is dependent on 
which face of the terrain that the corner of the rectangle lies in. It turns out that each possible comer 
placement can be handled individually without affecting the asymptotics. The n vertices impose an 
n X n grid on the domain, yielding O(n^) grid cells in which a comer may he. Because the terrain is 
a planar map, there are 0(n) edges as well, and each edge can cross at most 0(n) grid cells. Since 
no two edges can cross, this generates at most O(n^) new regions inside of all O(n^) grid cells. 
Since each rectangle is determined by the placement of two opposite corners the total complexity is 
not affected and is still O(n^). We summarize in the following lemma. 



reached. 



□ 



arg max 




M{p) — B{p) dp and arg max 
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Lemma A. 3. Consider a model where a system of A quadratic equations can be solved in 0{1) 
time. Then let h : M.'^ ^ be a piecewise-linear function with n vertices. 



can be solved in 0{n^) time. 



A. 3 Equations 

For a piecewise-linear terrain : — >^ R we wish to evaluate D{h,R) = Jjih{p) dp where 
R is some rectangle over the domain of h. Within R, the value of h is described by a set of 

triangles Tr = {ti, . . . , t^}. Let R be described by its four boundaries. Let xi and X2 describe the 
left and right boundaries, respectively, and let yi and 2/2 describe the top and bottom boundaries, 
respectively. Now 



Assume that we have computed the integral from xi up to x{v) the x-coordinate of a vertex v 
in the piecewise-linear terrain. To extend the integral up to the x{u) where u is the next vertex to 
the right. We need to consider all of the triangles that exist between x{v) and x{u). Label this set 
{ti,. . . , tk} where U is below tj in the y-coordinate sense for i < j. Note that no triangle can begin 
or end in this range and this order must be preserved. We also consider the intersection between the 
edge of the triangulation and R a vertex. Let the slope within a triangle ti be described 



arg max / 

RG0i2 Jr 





hi{x,y) = aiX + f3iy + ^i 



(A.l) 



and describe the edge of the triangulation that divides U and tj+i as 



1% — COiX -\- Ki. 



(A.2) 



Now the integral from x{v) to x{u) is described 




h{x, y)dydx 
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= / hi{x,y)dy + / hi{x,y)dy + / hk{x,y)dy 

Jx{v) Jyi i=2 'i-i(^) Jlk-iix) 



fx(u) 
v{v) 
rx{u) 
J x{v) 

x{v) 

^x(u) 
J x{v) 

J x{v) 



VI 

" ruJiX + Ki 

i-2 Jo 



dx 



{aix + Piy + ji)dy 

k—l fUliX+Ki 



{aiX + f3iy + ji)dy+ 



dx 



/c!lix+.._iK^ + + ^k)dy 
{aixyi + hPiivi)'^ + liVi) - {aix{ujix + ki) + + Kif + 'yi{uJix + ki))+ 



dx 



Eft— 
i=2 



k-1 f (a^x{u}i_ix + + + + 7i(a)j_i.T + Kj_i)) 



Q;ja;(Wi.T + K^) + ^(3i{L0iX + KiY + ^i{uJiX + Hi)) 



{akx{ujk-.ix + Kk-i) + i/3fe(a;fc_ix + k^-i)^ + ^ki^^k-ix + Kfc-i)) - {akxy^ + \l3k{y2Ylky2) 
{aixyi + \0iyl + 7iyi)+ 

\ Ei=i - Q^^; + (A+i - A)K^a;2 + 2KiWia; + nj) + (7,+! - 7i)- dx 

_ {0ikXy2 + iPkvi + lky2) 



\aiyix'^ + {\l3iyl + 7ijyi)a;+ 

k E-Ji -0i^ + 2K,uJi{pi+i - pi))x^ + " Pi))^'' + (7i+i - 7i + - 

lo:ky2x'^ - {,\(3kyl + iky2)x 

i(aiyi - aky2){x{uf - xivf) + {\Piyl - \Pkyl + liVi - lky2){x{u) - x{v))+ 



\{u^^,{(i,+r-(3,)){x{uf-x{vf)+ 
(7i+i - 7i + - Pi)){x{u) - x{v)) 



The point of all of this computation is that disc(/i, [xi, 2:2] x [yi, y2]) = /^^^ /^J^^^ ^(a^, y) dydx is 
a third order polynomial in X2, the x-position of the right endpoint, given that the set of triangles 
defining h is fixed. Thus to solve 



argmindisc(/i, [xi,X2] x [2/1,2/2]) 



X2 



requires finding where g|^disc(/i, [xi,X2] x [2/1,2/2]) = 0. If it is outside the range [x{v),x{u)], 
the two vertices bounding X2, then it is minimized at one of the two end points. It should be noted 
that by symmetry, the minimum of xi and X2 are independent, given 2/1 and 2/2, but that both are 
dependent on 2/1 and 2/2- Also, by symmetry, the previous statements can swap every x\ and X2 with 
2/1 and 2/2, respectively. 

Solving for argminj^.^ a,2]x[yi,j/2] disc(/i, [xi, X2] x [2/1,2/2]) requires solving 4 quadratic equations 
of 4 variables. This system is of the form 



1 To "j/i ] r a;i 1 To '^j/i ^^2/2 ] r ] r t^i 

^ I3y^ X2 ^ Ky^ Ky^ x\ ^ 7^^ 

a^^ a^.^ 2/1 ^^xi '^xi y\ -^y^ 

. J I Pxi Px2 J L 2/2 J I I^Xl 1^X2 J L 2/2 J L 72/2 

which in general has no closed-form solution. 
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B ^-approximations for a Normal Distribution 



A normal distribution, often referred to as a Gaussian distribution, is a family of continuous distri- 
butions often used to model error. The central limit theorem highlights its power by showing that 
the sum of independent and identically distributed distributions with bounded variance converge 
to a normal distribution as the set size grows. A normal distribution ^^(m, a^) is defined by two 
parameters, a mean m marking the center of the distribution and a variance scaling the spread of 
the distribution. Specifically, we can analyze a random variable X distributed according to a normal 
distribution, denoted X ~ J\f(m, u^). We then say that 

1 (x-rnf- 
(TV ZTT 

describes the probabiUty that a point X drawn randomly from a normal distribution ?\r(m, cr^) is 
located at x G M. Since it is a distribution, then J^^^ Vm,a'^ (^) = 1- A standard normal distribution 
!N(0, 1) has mean and variance 1. As the variance changes, the normal distribution is stretched 
symmetrically and proportional to a so the integral is still 1. The inflection points of the curve 
describing the height of the distribution are at the points m — a and m + a. 

A multivariate normal distribution is a higher dimensional extension to the normal distribution. 
A d-dimensional random variable X = [Xi, . . . , Xd]'^ is drawn from a multivariate normal distri- 
bution if for every linear combination Y = aiXi + . . .^a^Xd (defined by any set of d scalar values 
tti) is normally distributed. Thus for a cZ-dimensional random variable X defined over the domain 
W^, any projection of X to a one dimensional subspace of is normally distributed. 

We now discuss how to create e-approximations for (D, 3?2 x IK) where D is a multivariate normal 
distribution with domain R^. Extensions to higher dimensions will follow easily. We primarily 
follow the techniques outUned in Section 5.2 for a smooth terrain with properties z^, dxi, and Ad. 
What remains is to approximate D with another domain D' such that D has better bounds on the 
quantities z^, and di)'- The approximation will obey Lemma 5.1, replacing P with D' by just 
truncating the domain of D. 

The cumulative distribution function $^^2(0;) for a normal distribution (p^^a'^ describes the 
probabiUty that a random variable X ~ 3sr(m, cr^) is less than or equal to x. We can write 

where erf(x) = e~*^ dt. W.l.o.g. we can set m = 0. We want to find the value of x such 

that $0,0-2 (x) > 1 — e/4 so that if we truncate the domain of <y9o,<T2 (x) which is being approximated, 
the result will still be within e/2 of the original domain. We can bound erf(a;) with the following 
inequaUty which is very close to equality as x becomes large: 

l-erf(x) < 

x-v/tt 

For a a < 1/ {e^y^^) then 



1 — erf(a;) < a when x > y — ln(a\/7r)- 
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We say that the tail of ^I'o.cr^ is sufficiently small when 1 — $0,0-2 = ^ — ierf(x/((7\/2)) < e/4. 
Thus letting a = e/4, this is satisfied when 

x>a\j2hi{l/{e^f^)). 

Thus 



/„ 



m+CTiy21n(l/(£-yV2)) 

'y^'m,o-2 (a^) dx>l-e/2 



and bounding the domain of the normal distribution ^2 to 



m 



- (7-^2 ln(l/(£v/^)), + (7-^2 ln(l/(£v/^)) 



will approximate the distribution within e/2. 

For a multivariate normal distribution, we truncate in the directions of the x- and y-axis accord- 
ing to the variance in each direction. Letting T)' be the normal distribution with this truncated 
domain, then the diameter of the space is dxi' = 0(<7max-\/log(l/e)), where C7max is the maximum 
variance over all directions. (cTmin is the minimum variance.) In d-dimensions, the diameter is 

dq), = 0((Tmax\/c?log(l/£)). 

The lower bound z^, is now on the boundary of the truncation. In one dimension, the value at 
the boundary is 



( I , \ 1 -f(Tx/21n(l/(£^/W2))l /(2o2) £ 



For a 2-variate normal distribution the lower bound occurs at the corner of the rectangular boundary 
where in the 1-variate normal distribution that passes through that point and m = the value of 



-\/2(T Y 2 ln(l/ (£-^7r/2)). Thus the value at the lowest point is 

^D' = ^o,ct2 (^v2c7y21n(l/(£V7r/2))J = ^ V ) 



2c7c 

where is the variance in the direction of the comer. In the d-variate case, the lowest point is 

We calculate a bound for A©/ by examining the largest second derivative along the 1 -dimensional 
normal distribution and along an ellipse defined by the minimum and maximum variance. In the 
first case we can write 

dx^ = ^"''^^ 
which is maximized at x = \/3(7. And 



22 



For a bivariate normal distribution the largest eigenvalue of the Hessian of the extension of (pg ^2 
is similarly not large in the tail. Thus our choice of e does not effect this value. 

Hence, we can write {\d' d"^/ / Zx>') — Oi'p ^^S 7) for constant a. And, using Theorem 5.2, we 
can state the following theorem. 

Theorem B.l. For a 2-variate normal distribution with constant variance, we can deterministi- 
cally create an e-approximation of the range space {(p, x ^) of size O log^ ^ polylog(log ^)). 
This can be improved to a set of size 0{- log^ - polylog{\og -)) in time 0{\ polylog{log -)). 
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